hydrant <-
violation %>%
filter(violation %in% c("FIRE HYDRANT"), !is.na(geo_nyc_address))
pb <- progress_bar$new(total = nrow(hydrant))
get_coord <- function(url) {
pb$tick()
json_output <- fromJSON(url(url))$features
coord <- json_output$geometry[2]
borough <- json_output$properties$borough
out_df <-
tibble(
coordinates = coord$coord,
mapped_borough = borough
)
return(out_df)
}
hydrant_lat_long <-
hydrant %>%
mutate(
new_url = paste0("https://geosearch.planninglabs.nyc/v1/search?text=", geo_nyc_address, "&size=25"),
coord = map(new_url, get_coord)
) %>%
unnest(coord) %>%
filter(mapped_borough == borough)
hydrant_lat_long <-
hydrant_lat_long %>%
group_by(summons_number, geo_nyc_address) %>%
slice(1)
for (i in 1:nrow(hydrant_lat_long)) {
hydrant_lat_long$lat[i] <- hydrant_lat_long$coordinates[[i]][2]
hydrant_lat_long$long[i] <- hydrant_lat_long$coordinates[[i]][1]
}
write_csv(hydrant_lat_long, "hydrant_lat_long.csv")The previous chunk takes approximately 7 hours to run. We have saved the output from running this code in hydrant_lat_long.csv and the following code chunck should be run instead to load the data.
We will also load a dataset containing locations of all fire hydrants in NYC, from a file named Hydrants.csv.
hydrant_lat_long <- read_csv("hydrant_lat_long.csv") %>% select(-coordinates)## Rows: 71781 Columns: 34## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (17): plate, state, license_type, min, violation, precinct, borough, is...
## dbl (14): summons_number, hour, fine_amount, penalty_amount, interest_amoun...
## lgl (2): judgment_entry_date, coordinates
## date (1): issue_date##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.hydrant_actual <- read_csv("Hydrants.csv")## Rows: 109526 Columns: 9## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (2): the_geom, UNITID
## dbl (7): OBJECTID, BORO, POINT_X, POINT_Y, CB, LATITUDE, LONGITUDE##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.hydrant_to_plot <- hydrant_actual %>%
mutate(borough = "Hydrant") %>%
rename(
lat = LATITUDE,
long = LONGITUDE
) %>%
bind_rows(hydrant_lat_long)
hydrant_plot <-
hydrant_to_plot %>%
plot_ly(
lat = ~lat,
lon = ~long,
type = "scattermapbox",
mode = "markers",
alpha = 0.2,
color = ~borough)
hydrant_plot %>% layout(
mapbox = list(
style = 'carto-positron',
zoom = 9,
center = list(lon = -73.9, lat = 40.7))) We want to find the hydrants that are most commonly ticketed. As cars (or, drivers) are ticketed if they park within 15 feet of a hydrant, we will attempt to do this by constructing hitboxes of approximately 15 feet in radius around each hydrant and see how many violation points are within those circles.
Circles are a little difficult to work with. Let’s use squares instead, such that they are about 15 feet from the center of the square to any edge – this distance:
We used the following website to determine the distance each unit change in latitude/longitude covers in NYC. [-Link] (https://www.movable-type.co.uk/scripts/latlong.html)
If we go from the coordinate (40.7, -73.9) to (40.7, -74.0), we will have traveled 7.587 kilometers, or approximately 27657.48 feet. This indicates that a change of 0.000054 degrees longitude is approximately 15 feet.
If we go from the coordinate (40.7, -73.9) to (40.8, -73.9), we will have traveled 11.12 kilometers, or approximately 36482.94 feet. This indicates that a change of 0.000041 degrees latitude is approximately 15 feet.
To make things a little bit simpler, we will proceed as if one degree change in latitude covers the same amount of distance as one degree change in longitude. We will use the average of these two coordinate-distance approximations in our constructions of the hitboxes. We will just say that 15 feet can be represented by a change of 0.000045 degrees longitude or latitude.
# store "radius" of square (in degrees)
r = 0.00005
d = sqrt(2)*r
# need 5 points per square. points are: (x-d, y+d), (x+d, y+d), (x+d, y-d), (x-d, y-d), (x-d, y+d), (x-d, y+d)
working_fh <- hydrant_actual
square_coord <- function(xpt, ypt, dist = d){
x <- c(xpt - dist, xpt + dist, xpt + dist, xpt - dist, xpt - dist)
y <- c(ypt + dist, ypt + dist, ypt - dist, ypt - dist, ypt + dist)
out_mat <- cbind(x, y)
return(out_mat)
}
to_poly <- function(polymat, id){
poly <- Polygons(list(Polygon(polymat)), ID = id)
return(poly)
}
first_step <- working_fh %>%
mutate(
sq = map2(LATITUDE, LONGITUDE, square_coord),
polys = map2(sq, UNITID, to_poly)
) %>%
select(polys) %>%
pull(polys)
squares <- SpatialPolygons(first_step)
# now, get points ready
x = hydrant_lat_long %>% pull(lat)
y = hydrant_lat_long %>% pull(long)
xy = cbind(x,y)
dimnames(xy)[[1]] = hydrant_lat_long %>% pull(summons_number)
pts = SpatialPoints(xy)
# check if drawn squares have points in them
hits <- over(squares, pts)
hit_hydrants <- tibble(id = names(hits), hits = hits) %>% filter(!is.na(hits))
isolated_hydrants <- inner_join(hydrant_actual, hit_hydrants, by = c("UNITID" = "id"))isolated_hydrants_plot <-
isolated_hydrants %>%
mutate(borough = "Hydrant") %>%
rename(
lat = LATITUDE,
long = LONGITUDE
)
hydrant_plot <-
hydrant_lat_long %>%
plot_ly(
lat = ~lat,
lon = ~long,
type = "scattermapbox",
mode = "markers",
alpha = 0.1,
color = ~borough,
colors = "viridis")
plot <- hydrant_plot %>% add_trace(
data = isolated_hydrants_plot,
lat = ~lat,
lon = ~long,
type = "scattermapbox",
mode = "markers",
marker = list(
color = "red",
line = list(
color = 'white',
width = 2
),
size = 5
),
alpha = 1
)
plot %>% layout(
mapbox = list(
style = 'carto-positron',
zoom = 9,
center = list(lon = -73.9, lat = 40.7)))